%
% Complete documentation on the extended LaTeX markup used for Insight
% documentation is available in ``Documenting Insight'', which is part
% of the standard documentation for Insight.  It may be found online
% at:
%
%     http://www.itk.org/

\documentclass{InsightArticle}

\usepackage[dvips]{graphicx}
\usepackage{subfigure}

\usepackage{amsmath}
\DeclareMathOperator{\sgn}{sgn}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  hyperref should be the last package to be loaded.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[dvips,
bookmarks,
bookmarksopen,
backref,
colorlinks,linkcolor={blue},citecolor={blue},urlcolor={blue},
]{hyperref}


%  This is a template for Papers to the Insight Journal.
%  It is comparable to a technical report format.

% The title should be descriptive enough for people to be able to find
% the relevant document.
\title{Higher Order Accurate Derivative and Gradient Calculation in ITK}

%
% NOTE: This is the last number of the "handle" URL that
% The Insight Journal assigns to your paper as part of the
% submission process. Please replace the number "1338" with
% the actual handle number that you get assigned.
%
\newcommand{\IJhandlerIDnumber}{3231}

% Increment the release number whenever significant changes are made.
% The author and/or editor can define 'significant' however they like.
\release{1.0.0}

% At minimum, give your name and an email address.  You can include a
% snail-mail address if you like.
%\author{Galileo Galilei$^{1}$, Giordano Bruno$^{2}$ and Anthony Leeuwenhoek$^{3}$}
%\authoraddress{$^{1}$Pisa University, Tower drive\\
               %$^{2}$Rome University, Inquisition street\\
               %$^{3}$Netherlands Pragmatic University, Port street}
\author{Matthew McCormick$^{1}$}
\authoraddress{$^{1}$University of Wisconsin-Madison, matt@mmmccormick.com}

\begin{document}

%
% Add hyperlink to the web location and license of the paper.
% The argument of this command is the handler identifier given
% by the Insight Journal to this paper.
%
\IJhandlefooter{\IJhandlerIDnumber}


\ifpdf
\else
   %
   % Commands for including Graphics when using latex
   %
   \DeclareGraphicsExtensions{.eps,.jpg,.gif,.tiff,.bmp,.png}
   \DeclareGraphicsRule{.jpg}{eps}{.jpg.bb}{`convert #1 eps:-}
   \DeclareGraphicsRule{.gif}{eps}{.gif.bb}{`convert #1 eps:-}
   \DeclareGraphicsRule{.tiff}{eps}{.tiff.bb}{`convert #1 eps:-}
   \DeclareGraphicsRule{.bmp}{eps}{.bmp.bb}{`convert #1 eps:-}
   \DeclareGraphicsRule{.png}{eps}{.png.bb}{`convert #1 eps:-}
\fi


\maketitle


\ifhtml
\chapter*{Front Matter\label{front}}
\fi


% The abstract should be a paragraph or two long, and describe the scope of the
% document.
\begin{abstract}
\noindent In this article we describe higher order accurate derivative and gradient
image filters for the InsightToolkit.  These filters are central difference-based
numerical derivative approximations that account for additional Taylor
series terms and are based on the expressions given by Khan and Ohba.
\end{abstract}

\IJhandlenote{\IJhandlerIDnumber}

\tableofcontents

\section{Introduction}

A common way to approximate the first derivative of a function $f_k$ at sample offset $k$ using
finite differences is the central difference method.
\begin{equation}
  f^1_0  \approx \frac{f_1 - f_{-1}} { 2 h }
\end{equation}

Where $h$ is the sampling period.

This expression comes from a Taylor series expansion of the component terms
\begin{equation}
  f_1 = f_0 + h f^1_0 + \frac{h^2}{2!}f^2_0 + \cdots + \frac{h^n}{n!} +
  \mathcal{O}(h^{n+1})
\end{equation}

Where $\mathcal{O}(h^{n+1})$ indicates the series has been truncated after $n+1$
terms.

We also have
\begin{equation}
  f_{-1} = f_0 - h f^1_0 + \mathcal{O}(h^{2})
\end{equation}

Then we see
\begin{equation}
  f^1_0 = \frac{f_1 - f_{-1}} { 2 h } + \mathcal{O}(h^{2})
\end{equation}

Higher order accurate\footnote{Here we use the terminology \textit{order of
accuracy} to refer to the number of terms used in the Taylor series
approximation and \textit{order derivative} to refer to the degree of the
derivative.} approximations can be made by using additional samples.
For instance, a central difference approximation to the first derivative that uses a
five point kernel is
\begin{equation}
  f^1_0 = \frac{f_{-2} - 8 f_{-1} + 8 f_1 - f_2}{ 12 h } + \mathcal{O}(h^4)
\end{equation}

Khan and Ohba derived closed form expressions for higher order accurate
approximations for an arbitrary $p$th order derivative\cite{Khan1999,Khan2003}.

For example, the tap-coefficients, $d$, for a first order derivative 
\begin{eqnarray}
d_0 & = & 0 \\
d_k & = & (-1)^{k+1} \frac{N!^2}{k(N-k)!(N+k)!}, \; k = \pm 1, \pm 2, \dots, \pm
N
\end{eqnarray}

In this paper we present image derivative and gradient filters for the
InsightToolkit (ITK)\cite{itk} that implement these expressions.

\section{Implementation}

The interfaces of \code{itk::HigherOrderAccurateDerivativeImageFilter} and
\code{itk::HigherOrderAccurateGradientImageFilter} are similar to
\doxygen{DerivativeImageFilter} and \doxygen{GradientImageFilter}.  However,
these filters have an additional method, \code{SetOrderOfAccuracy()}.  The
approximation will be accurate in Taylor series terms to twice the
\code{OrderOfAccuracy}.  They are fast filters implemented with
\doxygen{NeighborhoodIterator}'s.  Note that the radius of the
\doxygen{NeighborhoodOperator} will be equal to the \code{OrderOfAccuracy}.
Each order derivative is implemented separately to use recursion during
calculation of the factorials\cite{Khan1999}.

It should be noted that the magnitude of this finite impulse response filter's components
falls off dramatically at the ends of the filter with increasing order accuracy.
Increasing the order of accuracy beyond four may not be beneficial.

\section{An Example -- Foot MRI Derivative and Gradient Magnitude}

\begin{figure}
  \begin{center}
    \includegraphics[width=5cm]{Comparisons/foot.png}
  \end{center}
  \caption{Foot magnetic resonance image (MRI) from which the derivative images
  in Figure \ref{fig:derivative} and the gradient images in Figure \ref{fig:gradient} are
  generated.  This image was taken from the ITK \textit{Testing/Data/Input}
  directory.}
  \label{fig:intensity}
\end{figure}

\begin{figure}
  \begin{center}
    \subfigure[DerivativeImageFilter.]{\label{fig:deriv-a}\includegraphics[width=5cm]{Comparisons/derivative_image_filter.png}}
    \subfigure[OrderOfAccuracy = 1.]{\label{fig:deriv-b}\includegraphics[width=5cm]{Comparisons/derivative_accuracy1.png}}
    \subfigure[OrderOfAccuracy = 2.]{\label{fig:deriv-c}\includegraphics[width=5cm]{Comparisons/derivative_accuracy2.png}}
    \subfigure[OrderOfAccuracy = 3.]{\label{fig:deriv-d}\includegraphics[width=5cm]{Comparisons/derivative_accuracy3.png}}
    \subfigure[OrderOfAccuracy = 4.]{\label{fig:deriv-e}\includegraphics[width=5cm]{Comparisons/derivative_accuracy4.png}}
    \subfigure[OrderOfAccuracy = 5.]{\label{fig:deriv-f}\includegraphics[width=5cm]{Comparisons/derivative_accuracy5.png}}
  \end{center}
  \caption[MRI foot derivative.]{Derivative along Direction 0 in the foot MRI
  image.
  The brightness in all the derivative images ranges from
  -110, black, to 110, white.  Figure \ref{fig:deriv-a} was created with the
  \code{DerivativeImageFilter} while Figure \ref{fig:deriv-b} to
  \ref{fig:deriv-f} were created with the
  \code{itk::HigherOrderAccurateDerivativeImageFilter} with the
  \code{OrderOfAccuracy} set from 1 to 5.}
  \label{fig:derivative}
\end{figure}

\begin{figure}
  \begin{center}
    \subfigure[GradientImageFilter.]{\label{fig:grad-a}\includegraphics[width=5cm]{Comparisons/gradient_image_filter.png}}
    \subfigure[OrderOfAccuracy = 5.]{\label{fig:grad-b}\includegraphics[width=5cm]{Comparisons/gradient_accuracy5.png}}
    \subfigure[DoG.]{\label{fig:grad-c}\includegraphics[width=5cm]{Comparisons/dog.png}}
    \subfigure[Recursive Gaussian.]{\label{fig:grad-d}\includegraphics[width=5cm]{Comparisons/recursive_gaussian.png}}
  \end{center}
  \caption[Gradient magnitude.]{Gradient magnitude of Figure \ref{fig:intensity}
  with different filters.  The brightness in all the gradient images ranges from
  0.0, black, to 128.0, white.  Figure \ref{fig:grad-a} is the
  \doxygen{GradientImageFilter}.  Figure \ref{fig:grad-b} is the
  \code{itk::HigherOrderAccurateGradientImageFilter} with order of accuracy set
  to 5.  Figure \ref{fig:grad-c} is the
  \doxygen{DifferenceOfGaussiansGradientImageFilter} with a width of 1.0.
  Figure \ref{fig:grad-d} is the \doxygen{GradientRecursiveGaussianImageFilter} with
  a sigma of 1.0.  The spacing in the image is isotropic with a value of 1.0.}
  \label{fig:gradient}
\end{figure}

Figure \ref{fig:derivative} and Figure \ref{fig:gradient} demonstrate the behavior of the filter on an MRI foot
image.  Note that if the \code{OrderOfAccuracy} is set to one,
\code{itk::HigherOrderAccurateDerivativeImageFilter} is equivalent to
\doxygen{DerivativeImageFilter} and
\code{itk::HigherOrderAccurateGradientImageFilter} is equivalent to
\doxygen{GradientImageFilter}.  The source code to generate
these images can be found in the \code{Testing} subdirectory.

In general, we see the higher order accuracy results in lower noise, higher
contrast and higher dynamic range results.  Unlike the difference of Gaussians
or recursive Gaussian filters, though, it does not remove high frequency
content, which is often primarily noise.  Choice of the most appropriate filter
depends on the application.

\section{Future Work}

At this time, only the first order derivative is implemented.  In the future,
higher order derivatives should be implemented.

\section{Acknowledgments}

The authors would like to thank the NIH for their funding with grants
T90DK070079 and R90DK071515.

\bibliographystyle{itk}
\bibliography{ij_higher_order_accurate_gradient}


\end{document}
